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Abstract. The Accretion- E jectio n Instability (AEI) de- 
scribed by Tagger & Pellat (jl99£| , hereafter TP99) is ex- 
plored numerically using a global 2d model of the inner 
region of a magnetised accretion disk. The disk is initially 
currentless but threaded by a vertical magnetic field cre- 
ated by external currents, and frozen in the flow. In agree- 
ment with the theory a spiral instability, similar in many 
ways to those observed in self-gravitating disks, develops 
when the magnetic field is, within a factor of a few, at 
equipartition with the disk thermal pressure. Perturba- 
tions in the flow build up currents and create a perturbed 
magnetic field within the disk. The present non-linear sim- 
ulations give good evidence that such an instability can 
occur in the inner region of accretion disks, and gener- 
ate accretion of gas and vertical magnetic flux toward the 
central object, if the equilibrium radial profiles of density 
and magnetic flux exceed a critical threshold. 
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1. Introduction 

Magnetic fields in accretion disks have been seen over the 
last decade to strongly effect the structural and dynamical 
properties of accretion disks as compared to the well de- 
fined standard steady-state, thin disk models (Shakura & 
Sunyaev |1973[ , Pringle [1981[ Frank, King & Raine |1992| ). 
Since the early 1990s, numerical simulations of accretion 
disks have increasingly shown huge deviations from the 
smoothly varying, symmetrical (azimuthally and verti- 
cally), time-independent accretion disk models. 

Although the magnetic field has been assumed to be 
of importance in generating instabilities within the disk 
(leading to turbulent transport), this has essentially been 
modeled as turbulent viscosity in the classical prescription 
v = ac s h where v is the viscosity, c s the sound speed, h the 
disk height and a a dimensionless parameter. This then 
allows for the outward radial advection of angular momen- 
tum and inward radial transport of matter. However, this 
alone is insufficient to explain the time- variability of accre- 



tion disks such as relatively long timescale (hours to days) 
dwarf novae ou tburst of cataclysmic variables (e.g. Can- 
nizzo & Mattei |l992|) or mor e rap id (of the order of sec- 
onds) QPO s (e.g. Swank et al \L99Tj Munoef al. § Sobczak 



et al. 2000 ). The time variability of these systems has al 
ways been shown to be associated with the dynamics of 
the accretion disk itself - a feature that is not included 
within standard hydrodynamical treatments. 

Since the discovery of the importance of the magneto- 
rotatio nal (M RI) or Balbus- Hawlcy (BH ) instability (Ve- 
likhov 1959| , Chandrasekhar 1960 1961 ) in the context of 
accretion disks by Balbus & Hawley ( 1991), and its numer- 
ical demonstration by Hawley & Balbus ( 1991 ), a number 
of local 3d simulations have been performed (e.g. Bran- 
denburg et al. 1995, Hawley et al. 1995, Stone et al. 1996, 



Hawley et al. [1996| , Hawley |2000b| ). These all show that 
turbulent transport is indeed generated through a local in- 
stability. The resulting transport, if it is measured in terms 
of the a prescription, is in the range expected from ob- 
servations; however the resulting a parameter is strongly 
varying, meaning that the turbulence does not obey the 
simple prescriptions of the a model. Other features, such 



as time variability (Hawley & Krolik 2000), a vertically 
varying a (Brandenburg et al. 1996b) and vertical struc- 
tural asymmetries (Caunt |l998| , Miller & Stone 2000), im- 



ply that the magnetic field strongly affects the overall dy- 
namics of an accretion disk (and in particular its vertical 
structure and its connection to the corona). 

Jets, associated with accretion disks in Active Galac- 
tic Nuclei, X-ray binaries and Young Stellar Objects are 
certainly magnetised (Blandford & Payne 1982) and the 
magnetic field at the disk surface is required to drive the 
matter across the slow magnetosonic point, after which 
centrifugal forces take over (Lovelace et al. 1991). The 



formation of jets is also observed to be connected to the 
accretion within the disk. In fact one the main attrac- 
tions of MHD jet models is that they show that the jet 
is a very efficient way to extract angular momentum from 
the disk, in sharp contrast with a-disk models where the 
angular momentum is evacuated outward radially by vis- 
cous stresses. A number of instabilities other than the MRI 



have been proposed for accretion disks (e.g. Parker 1979 
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Papaloizou & Pringle 1984, Spruit et al. 1995); however 
the recent discovery of the accretion-ejection instability 
(AEI) (Tagger et al. |1990| , TP99) provides a possible link 
between radial and vertical angular momentum transport. 

As shown by TP99, the instability occurs close to the 
inner edge of an accretion disk for which the quantity 
0£/£>o increases radially where Vt is the angular velocity, 
£ the surface density and Bq the magnetic field for which 
the plasma beta (the ratio of gas to magnetic pressures) 
is of a moderate value (around unity). The instability ap- 
pears as a low azimuthal wavenumber spiral wave, made 
unstable through the interaction with a Rossby vortex (as- 
sociated with the gradient of vorticity) which it generates 
at its corotation radius. This also causes the emission of 
Alfven waves from corotation vertically along magnetic 
field lines, as their footpoints are twisted. Hence, in the 
context of jet formation, the AEI naturally contains the 
highly desirable property of linking vertical propagation 
directly to the disk dynamics. Also, unlike the local MPJ 
which is quenched for moderate fields, the AEI exists in 
situations for which the field pressure is comparable to 
the gas pressure, as may occur in the inner region of an 
accretion disk as magnetic field is advected inwards with 
matter and accretes towards the central object. 

With the view of the low-frequency QPOs observed in 
X-ray binaries, in particular those hosting a black hole, the 
instability is also a good theoretical candidate (Varniere 
et al., 200C; Rodriguez et al, |2000p . One can thus imagine 
a scenario in which, under initial low field strength situ- 
ations the MRI initially acts, advecting the field radially 
inward by turbulent transport. If at least part of the ver- 
tical magnetic flux is advected with the gas the field in 
the inner disk region builds up, until a certain point is 
reached where the higher magnetic pressure quenches the 
MRI an d the AEI sets in. This 'magnetic flood scenario' 
(Tagger 1999 ) might explain the ~ 30 minutes cycle s of 
the micro-quasar GRS 1915+105 (Mirabel et al. |l99Sj ). 

The inner region of accretion disks is increasingly be- 
coming the focus of 3d numerical experiments. As com- 
putational power increases, global models (still limited by 
resolution and radial extent) are now being produced as 
well as t he we ll-established local shearing-box models. Ar- 
mitage ( |l998|) describes a global model that omits vertical 
stratification and is limited in both vertical and radial ex- 
tent, but illustrates effectively th e globa l development of 
the MRI More recently Hawley (|2000aJ ) and Machida et 
al. (|2000|) describe simulations of accretion disks starting 
from toroidal configurations (i.e. a departure from Kep- 
lerian thin disks); at the time of preparation Hawley & 
Krolik ( 2000 ) and Hawley ( 2000b ) have published numer- 
ical simulations of full accretion disk models. 

We note here, and will discuss in our conclusions, that 
in these 3d simulations spiral waves are seen to occur close 
to the inner boundary of the disk as a result of magnetic 
stresses; this is very similar to the features expected from 



the AEI, and indeed very similar to the results presented 
here. 

In a very different approach, Stehle & Spruit 200C have 
presented a 2D (r, <fi) model of a disk. The disk is con- 
sidered as infinitely thin and embedded in vacuum. Their 
model is used to demonstrate the existence of the inter- 
change instability, predicted to occur in strongly magne- 
tised disks (Spruit & Taam 1990, Spruit et al. 1995), and 
associated with the radial gradient of the magnetic field. 
Unexpectedly, they find that the disk is also unstable for 
lower field strength than predicted. 

Our model is very similar to theirs, with numerical 
characteristics optimised according to our previous knowl- 
edge of the AEI. We also consider an infinitely thin disk 
threaded by a moderate vertical magnetic field: this is jus- 
tified by the properties of the AEI, which is essentially con- 
stant vertically across the disk. The disk is also embedded 
in vacuum but, in order to separate different physical ef- 
fects, we consider in the present paper only configurations 
where initially the equilibrium magnetic field is due to ex- 
ternal currents. As discussed below we use a logarithmic 
radial grid: this allows us to get a better resolution in the 
inner region of the grid, where it is necessary, and on the 
other hand to model a much larger radial extent of the disk 
and to get rid of unwanted effects of the boundary condi- 
tion at the outer radius of the simulation. We believe that 
such effects may explain the stronger instability observed 
in their work, compared with ours. This numerical setup 
(with, as discussed below, a magnetic potential given by a 
Poisson equation) allows us to use well-established meth- 
ods developed in the description of self-gravitating disks. 

In spite of these differences, our results arc quite sim- 
ilar to the ones they obtain at moderate field amplitude 
and we believe that, although the initial conditions dif- 
fer and the numerical methods are independent, similar 
physics are at the heart of the two simulations. 

This model allows for the development of the spiral 
wave and the vortex - at the heart of the instability. The 
paper is organised as follows: Section 2 describes the basic 
model including the fluid equations and implementation of 
the magnetic field, Section 3 provides details of the numer- 
ical methods used, in Section 4 we present our results and 
in Section 5 our conclusions. 



2. The model 

We present, in cylindrical coordinates, a global two- 
dimensional model of the inner region of an accretion disk 
around a central object. The model is 2D in the sense that 
we solve only for the horizontal (r, </>) components of the 
velocity, and that all perturbed quantities are assumed in- 
dependent of z within the disk. We assume that the disk is 
threaded by a poloidal magnetic field which is symmetric 
about the midplane, hence purely vertical at z = 0. As 
described later, assuming that the disk is of a finite ver- 
tical height and embedded in vacuum, the magnetic field 
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outside the disk (and in particular at its surfaces) can be 
derived from a magnetic potential, which in turn, is cal- 
culated from the distribution of magnetic field within the 
disk (this is analogous to calculating the gravitational po- 
tential in a self-gravitating disk). Currents within the disk 
can be derived from the jump in the horizontal component 
of the field across the disk. Hence magnetic stresses can 
be included in an otherwise purely hydrodynamic model. 

The variables we solve for are the radial and per- 
turbed azimuthal velocity, u r and u$, the surface density, 
E (which comes from the vertical integration of the den- 
sity, p) and the vertical component of the magnetic field, 
B z . The magnetic potential is solved separately as <& m . 
The "perturbed" component of the azimuthal velocity is 
defined as its deviation from keplerian rotation (so that it 
contains both the true perturbation and the equilibrium 
contribution of the pressure gradient). By substituting for 
the full azimuthal velocity t/^ = uk + u^, where uk is the 
Keplerian velocity, the axisymmetric component of U<f> can 
be enforced more accurately and is less prone to numerical 
roundoff errors (encountered from large azimuthal veloci- 
ties and gravitational forces) which can lead to anomalous 
radial motion. 

The instability can be successfully described using a 
2d model because the vertical scale height of the pertur- 
bations is greater than the disk thickness (TP99), unlike 
the Magneto-Rotational Instability (MRI) for which the 
vertical wavelength is critical, so that 3d local models are 
required to study the full nature of the instability {e.g. 
Brandenburg et al. 1995, Hawley et al. 1996, Stone et al. 



1996). Hence our simulations will not show the full tur- 



bulent transport of angular momentum expected within 
accretion disks, but will however allow us to isolate the 
effects of the AEI alone in their inner region, providing 
a first approach of the physics of the AEI, and valuable 
details for future 3d models. 

2.1. Fluid equations 

The model solves for the momentum, induction and con- 
tinuity equations integrated vertically across the disk, in 
the hypothesis that (as expected for the AEI) the vertical 
scale height of the perturbations is larger than the disk 
thickness. The momentum equations for r and (f> are given 

by 



du r 
~dt 



du^ 
dt 



du r U<f, du r u $ 2u c / ) uk 
dr r dej) r r 
1 dP B z AB r 
~E d7 + 4ttE ' 



dr 



(1) 



r c 
1 dP 
Er Ikb 



4ttE 



(2) 



where E is the surface density, = + and = 
y/GAl/r is the Keplerian velocity. Since pressure stresses 
are relatively unimportant, we assume for simplicity an 
isothermal equation of state so that the pressure gradient 
only changes as a result of advection of fluid density. The 
vertically integrated pressure, P, is given by 



P = c E, 



(3) 



where c s is the sound speed. The terms AB r and AB^ 
represent the currents associated with the jumps in radial 
and azimuthal magnetic field across the surfaces of the 
disk. Assuming symmetry conditions apply then AB r = 
B+ — B~ = 2B+ where the labels + and - note the values 
at the upper and lower disk surfaces, and similarly AB^ — 

The continuity equation, integrated vertically, gives 



9E 
dt 



ld_ 

r dr 



(Eru r ) 



(4) 



Similarly the vertical component of the induction equation 
takes the form: 



dB z 
dt 



ld_ 

r dr 



(B z ru r 



^(B z r U4> ). 



(5) 



The identical form of the continuity and induction 
equations results from flux freezing in the ideal MHD 
limit, i.e. from the fact that both surface density and ver- 
tical magnetic flux are conserved in each tube threading 
the disk. It will be used below in the calculation of the 
magnetic potential. 

2.2. Magnetic potential 

Equations |J and [| show that the horizontal advection of 
the vertical magnetic field has an identical form to the 
advection of density in the limits of a vertically averaged 
model. Another identical relationship can be derived from 
assuming that the disk lies in vacuum. With this condi- 
tion the absence of currents outside the disk allows one to 
describe the perturbed magnetic field outside the disk as 
the gradient of a magnetic potential, $m- Following the 
derivation in TP99, we describe the field outside the disk 



as: 



B = -sign(z)V$ 



(G) 



so that both B z and <&m are even in z, while B r and 
B,), are odd. In vacuum the divergence-free field condition 
gives 



V 2 $ju = 0. 

On the other hand symmetry gives: 

— $ M (z = 0+) = - — $ M (z = (T) = —B 1 
dz dz 



(7) 



(8) 
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where is the perturbed value of B z in the disk. Thus 
we can write throughout space: 



V 2 $ 



M 



-2B°5{z), 



(9) 



where S(z) is the Dirac function. The form of this is identi- 
cal to the standard Poisson equation for the gravitational 
potential of a self-gravitating disk: 



(10) 



Hence, from the similarity between the Poisson equations 
on one hand, and on the other hand the similarity be- 
tween the continuity equation (^) and the induction equa- 
tion (^) , we can solve for the magnetic potential using the 
classical Poisson kernel and more generally the methods 
devised in the simulation of self-gravitating disks (Binney 
& Tremaine 1987), as described in Section 3.3. The inte- 
gration of the potential from individual field lines within 
the disk yields an expression for the total magnetic poten- 
tial as: 



T2 p2tt 



2tt J ri J Q 



B^(r',eydr'd9 



\/ r 2 + : 



2rr' cos 9 



(11) 



This similarity explains that, as found by Tagger et al. 
( |1990D , the waves we will describe are very similar to the 
spiral waves of self-gravitating disks. An important differ- 
ence, however, is the negative sign in Equation meaning 
that the magnetic field plays the role of a negative self- 
gravity: its presence "stiffens" the plasma against com- 
pression, whereas self-gravity helps it. 

2.3. Initial Conditions 

For the models presented here we assume, in the context 
of the accretion disk of a black hole in an X-ray binary, 
that the central object has a mass of M = 10M Q and the 
disk extends from between r = 1, 000 km to r = 50, 000 km 
from the origin at r — with a constant aspect ratio of 
e = h/r = 0.1 where h is the disk thickness. The only crit- 
ical parameters, however, are e and the ratio (3 of thermal 
to magnetic pressure, so that similar results would apply 
in disks of similar aspect ratio but very different dimen- 
sions, such as the accretion disks of young stellar objects 
or Active Galactic Nuclei. 

The sound speed within the disk is calculated directly 
from the thin-disk approximation (e.g. Frank, King & 
Raine 1992) such that c s = fiflx where £Ik is the Keple- 
rian angular velocity; hence c s varies radially as r -1 / 2 . The 
disk starts from a steady-state velocity profile to which we 
add small scale random fluctuations with, unless otherwise 
stated, a maximum Mach number of 0.1 (hence the abso- 
lute magnitude of the initial fluctuations also decreases 
radially). The initial steady-state azimuthal velocity is set 
from the equilibrium condition, such that the centrifugal 
force counteracts radial forces acting on the fluid which 
include both gravity and gas pressure. Initially there are 



no currents within the disk and the magnetic field is of a 
purely external origin, being purely vertical with a radial 
gradient. This results in magnetic stresses which are ini- 
tially vanishing, but build up with the magnetic potential 
during the simulation. The exact profiles of B z , the ini- 
tial external field, and S, the surface density, are varied to 
examine the criterion of instability. Linear theory (TP99) 
lets us expect the disk to be unstable for initial profiles 
such that the quantity il£/£? 2 in crea ses radially. This cri- 
terion is investigated in Section |4.3| . The strength of the 
external field is also a variable parameter. The plasma 
beta, using the vertically integrated quantities, is defined 
as 



8ttP %ttY<c s Q k 



B 2 



B 2 



(12) 



taking into account standard thin disk approximations. 
We vary the field strength using different values of (3 de- 
fined at the inner edge of the disk. The instability occurs 
for moderately strong magnetic fields and is theoretically 
most active for (3 ~ 1. We therefore perform tests on f3 in 
the range 0.1 < < 10. 

2.4. Non-dimensionalisation 

We choose typical length, time and density scales to 
non-dimensionalise the quantities. We take [l]=10 8 cm = 
1000km, [t]=1000s and [p] = lO^gcuT 3 . These have 
the added advantage that velocities are returned as [v] = 
lkms -1 . The disk dimensions are then r = {1,50}, the 
surface density, we take to be S = 10[p][i] = lkgcm~ 2 
and the inner orbital period is then T r=1 = 1.7 x 10 _4 [£] = 
0.17 s. 



3. The numerical methods 

Equations to || are solved on a 2d cylindrical (r, 4>) grid. 
A conser vative scheme of the type described by Stone & 
Norman ( |l992[ ) is used to advance the fluid properties in 
time. The scheme uses a staggered mesh such that £ and 
B z are cell centered quantities and u r and are edge 
centered. This is advantageous in the fact that deriva- 
tives, required for a particular variable, are typically calcu- 
lated at the desired position and with increased accuracy. 
The grid itself is evenly spaced in azimuth and logarithmi- 
cally in radius. This has two advantages for this particular 
model: the region of interest (close to the inner bound- 
ary) has greater resolution and waves propagating to the 
outer boundary are damped via inherent numerical dissi- 
pation, which guarantees the outer boundary condition for 
these simulations (radiation condition, i.e. no information 
coming in from the outer boundary). This is particularly 
important here in reference to the related Papaloizou- 
Pringle instability (Papaloizou and Pringle, 1984), which 
we should recover in the limit of vanishing magnetic field 
(/3 — * 00). In this limit it is well known that a reflecting 
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outer boundary can result in a much more violent insta- 
bility (Narayan et at, 1987). We believe that this might 
explain the apparently more unstable numerical results of 



where va = B^er/8wT, is the Alfven velocity. The mag- 
netic contribution to the timestep is therefore of the order 
of 



Stehle and Spruit ( pOOOQ . 

Averaging of quantities onto the staggered mesh takes 
into account the non-uniform nature of the grid. The res- 
olution is typically taken to be 256x128 cells in r and <fi 
although lower resolution tests have been made for com- 
parison. The high radial resolution is used to capture the 
fine details of the instability. The azimuthal resolution is 
arguably excessive since the instability generates low az- 
imuthal wavenumber spiral waves. However we found in- 
accuracies, or even numerical instabilities, associated with 
the regularisation of the logarithmic singularity of the 
Poisson kernel described by Binney & Tremaine ( 1987 ) 
when the azimuthal and radial grid sizes are very different. 
Tests with modified smoothing schemes have been made 
with lower azimuthal resolution and indicate that the re- 
sults presented here are not dependent upon the scheme, 
or higher resolution. Further details of the Poisson kernel 
are given in Section 3.3. 



3.1. Evolution of the fluid 

The equations are solved in a conservative form such that 
angular momentum, mass and total magnetic flux are con- 
served to numerical precision throughout the simulation. 
The method follows that of Stone & Norman (1992) such 
that the terms are split into source and transport terms: 
the source terms being centrifugal, gravitational, pressure 
and magnetic terms in the momentum equation and the 
transport terms being the advection of matter. 

The source terms are calculated initially using stan- 
dard second order finite difference operators on the stag- 
gered mesh (modified for a logarithmically spaced grid) 
followed by the transport terms. The transport ter ms are 
calculated using a standard second order Van Leer ( 1979 ) 
upwind scheme. 

We also invoke the FARGO scheme (Masset |2000f) to 
allow for increased timesteps, typically limited by the un- 
derlying large velocity Keplerian motion. The Courant- 
Friedrichs-Lewy condition on the maximal timestep for 
numerical stability is thus reduced from dt £ dx/vg to 
dt £ dx/vg, where dx is the grid size and v g the group 
velocity of waves which can propagate in the simulation 
box. Typically v g will be of the order of the fast mag- 
netosonic speed, i.e. comparable with the sound speed if 
(3 ~ 1. The allowed timestep is thus increased by a fac- 
tor ve/vg ~ e , since vg is dominated by the keplerian 
velocity. 

In fact v g can be estimated more precisely: from the 
dispersion relation given by TP99, the maximum group 
velocity can be approximated as 



(13) 



dt = min ( Ax^r- 

V A 



taken over the whole of the numerical domain. 



(14) 



3.2. Boundary conditions 

The effect of boundary conditions on the behavior and 
development of the AEI was discussed in detail in TP99. 



3.2.1. Inner boundary 

Concerning the inner boundary, they showed that any re- 
flecting condition {e.g. u r = 0, ug = 0, etc.) results in the 
existence of modes with the same instability criterion, but 
different frequency. A detailed discussion of the physics of 
the inner boundary of the disk is impossible: indeed, in 
the example of the accretion disks of black holes in X-ray 
binaries, we still do not know what exists between the disk 
and the black hole horizon; one could consider an ADAF, 
or the magnetic structure threading the black hole, asso- 
ciated with the Blandford-Zjanek (1977) mechanism. In 
both cases the detailed physics of the inner disk boundary 



is unknown. Hawley and Krolik (2000) and Hawley( 2000b) 
have recently used a pseudo-newtonian gravitational po- 
tential to model the flow of the gas at the Last Stable Orbit 
near a black hole. This removes the need for an artificial 
boundary condition at the inner disk edge. However the 
low-frequency QPO which we consider as a manifestation 
of the AEI in X-ray binaries is observed mostly when the 
dis k stay s away from the Last S table Orbit (Varniere et 
al., 200C; Rodriguez et at, 2000 ). Here for simplicity we 



adopt the commonly used condition u r = 0. 

A disadvantage of this condition is that, since mag- 
netic flux is advected with the gas, it piles up in the inner 
few grid-zones. However, we find that overall, the effective 
change in the field strength is acceptable for our current 
models, where we find that the plasma (3 changes by 70% 
over the course of 200 orbits in the most extreme case. We 
feel that, although this is not the most desirable effect, it 
is acceptable for the purposes of this paper. Future mod- 
els will allow for the magnetic flux to spread in the inner 
"hole", between the disk and the central object. This will 
limit the build-up of the magnetic field (and its gradient) 
near the inner disk boundary. 

3.2.2. Outer boundary 

At the outer boundary, reflecting conditions are often used 
for simplicity but are non-physical, unless one can argue 
for the existence of a sharp outer disk radius. For instance 
in the case of an X-ray binary, a sharp outer boundary 
may exist, associated with the physics of tidal effects and 
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the gas flow from the companion. We have already men- 
tioned that, as shown by Narayan et al. (1987), a reflecting 
outer boundary may make the instability much more vio- 
lent. However in realistic conditions the ratio of the outer 
to inner disk radius is very large, and this effect should not 
act. Here the logarithmic radial grid allows us to extend 
the simulation grid to very large radius, and there the grid 
size also becomes large. On the other hand for modes of 
interest the dispersion relation (Tagger et at, 1990) shows 



that the radial wavelength becomes small at large radius. 
Numerically, when it becomes comparable with the grid 
size, the waves are heavily damped (this is a natural conse- 
quence of the use of low order derivatives and advection). 
This ensures that waves generated in the region of physi- 
cal interest, at small or moderate radius, will be damped 
before they reach the outer boundary, and are thus not re- 
flected. This effectively implements the correct boundary 
condition that no information should flow from infinity. A 
further advantage of the large radial interval is also that, 
for a double check, we can push the outer boundary far 
enough that waves cannot travel to the outer boundary 
within the limited duration of the simulations. 

3.3. Calculation of the magnetic potential 



As discussed in Section 2.2, the magnetic potential can be 
solved analogously to the gravitational potential in a self- 
gravitating disk. For speed of calculation w e use the FFT 
method discussed by Binney & Tremainc ( 1987 ). Rather 
than using the perturbed density, the perturbed vertical 
magnetic field is used as the source of the potential, and 
the Poisson kernel is modified to account for the change in 
coefficients (and sign) between Equations || and M. As in 
Binney & Tremaine, we perform a change of variables to 
reflect the logarithmic grid and simplify the calculation. 
It can thus be shown that after transforming the magnetic 
field and potential such that 



V M {u,<t>) = exp(u/2)$ M [r(u),<f>] 
b z (u,<t>) = exp(3u/2)Sf (r(u),<f> 



(15) 
(16) 



where Vm is the reduced potential, b z the reduced per- 
turbed magnetic field and u the logarithm of the radius, 
r, Equation [ll] can be rewritten as 

VmM) = — \ 



where Au and A<p are mesh sizes. However this form is 
only valid for small and comparable values of Au and A<f>. 
Here the physics of the AEI lets us expect (as confirmed 
by the results below) only low-m instabilities (where m is 
the azimuthal wavenumber), and for speed of computation 
we often wish to use much fewer grid cells in <f> than in it. 
In that case the approximation ( [l9| ) gives the unphysical 
property that JC(0,0) < K(Q, 1), which in turn generates 
numerical instabilities. In this case we use an alternate 
method, often used in self-gravitating disks, modeling the 
finite thickness of the disk (which does physically regu- 
larise the Poisson kernel). In order to do this we replace, 
in the computation of the kernel, the denominator \r — r'\ 
by yj\{r — r') 2 + h 2 . This makes the kernel regular at 
r = r'; it removes the need for a separate evaluation of 
K(0,0), which can be written as 

K(Au, A(f>) = 1 (20) 

^2[cosh(Au) - cos(Ac/>) + e 2 exp(Au)] 

where Au = u—u', Acf> = 4>—4>' and e is the aspect ratio of 
the disk. For the purposes of the magnetic potential alone, 
we have experimented with modifying e to effectively al- 
ter the smoothing length (however keeping the disk height 
identical) and find that it makes only a very slight differ- 
ence to the extent of the disk that is affected strongly by 
the instability. The overall development of the instability 
is unchanged. This was to be expected since the smoothing 
affects only short-range interactions, while the instability 
has large wavelengths except at large radii where we wish 
it to be damped by the limited resolution of the grid. For 
the purposes of this paper, in showing the general charac- 
teristics of the instability, we generally use the Binney & 
Tremaine kernel; however we include later an example of 
a low azimuthal resolution simulation using the modified 
kernel ©. 

The power of the FFT method of calculation is that 
the kernel (and its Fourier transform) need only be cal- 
culated once at the beginning of the simulation. At the 
beginning of each time step the Fourier transform of b z 
is calculated. The inverse Fourier transform of the con- 
volution of this and the Fourier transform of K is then 
performed to generate Vm which in turn gives $m- 

The radial and azimuthal magnetic field at the disk 
surface can then be calculated from 



K(u-u' 7 <j)-(f>')b z (u',(j)')du / d<j)' {17) B = -sign(z)V$ 



with the Poisson kernel defined as 



K{u-u',cj)-(j)') 



-y/2[cosh(w — u') — cos(0 — </>')] 



(18) 



This kernel has a logarithmic singularity where u — u' = 



and (f> — (/>' = 0. For this case Binney & Tremaine (1987) 
give an approximate smoothing: 



#(0,0) 



— — arcsinh ( -t—\ 
Ac/) \AuJ 



— — arcsinh ( -— - ] 
Au \A(j) J 



(19) 



m (21) 

yielding B r and B^. Their jump across the disk (coming 
from the sign of z in Equation [2l]) gives the current in 
the disk, and thus the Lorentz force in the source terms 
of Equations [j] and ||. 

The initial setup of the field is somewhat unphysical 
since it assumes that the field is initially totally of external 
origin, and we have no currents in the disk. The current 
density in the azimuthal direction is given by 

H = — —=°> ( 22 ) 



dz 



Or 
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Fig. 1. Comparison of the evolution of instability for different initial values of (3 = 8np/B 2 , the ratio of thermal to 
magnetic pressure. Values of /3 = 0.1, (3 = 0.5, f3 = 1.0, (3 = 5.0 and f3 = 10.0 are shown as dotted, dot-dashed, 
short-dashed, dot-dot-dot-dashed and long-dashed lines respectively. Time is given in units the orbital time at the 
inner radius. The plots show the evolution of (top) the random velocity, and (bottom) the average magnetic field. 



therefore, for a magnetic field which is purely vertical at 
the disk, and decreases radially one has 



(Mr 

dz 



< 0. 



(23) 



This implies that the field lines must bend inward towards 
the central object; in a physically realistic situation the 
opposite is usually considered, giving the familiar "hour- 
glass" shape commonly used in disk-jet models. However, 
for the purposes of this paper, we prefer to use this setup 
which allows us to identify the effect associated with the 
field gradient, and leaves out the action of field curvature 
and equilibrium currents. Preliminary linear calculations 
have shown us that an "hourglass" configuration is slightly 
more unstable than the one we use here. The code is per- 
fectly able to describe such initial configurations, with cur- 
rents both in the central hole and within the disk. These 
will be considered in future publications. 

4. Results 

Our main goal is to characterize the instability, its depen- 
dence on the properties of the disk, and its influence on 
the global evolution and accretion rate. We will thus first 
vary different parameters affecting the instability, and its 
strength. We will then discuss in selected cases the de- 
tailed properties of the instability, and finally evaluate its 
efficiency for accretion. 

4-1- Effect of the field strength 

Let us first consider the influence of the magnetic field 
pressure in the disk. For this we choose equilibrium mag- 
netic field and density profiles which, according to the 
linear theory, should give a rather strong growth rate: the 
instability criterion (TP99) is that 



d qe 

— In —J- > 0. 

or Bq 



(24) 



We thus take an equilibrium magnetic field varying as 



B (r) = Bo(rmin) — 
r 



5/4 



(25) 



which is the profile used by Stehle and Spruit ( |2000D . This 
should be unstable if the surface density profile is flatter 



than £ 



For this first series of runs we choose a 



somewhat unrealistic flat density profile, and vary the pa- 
rameter (3. We know that the magneto-rotational instabil- 
ity exists only for (3 > 1, whereas the AEI is most unstable 
for /? ~ 1. 



Fig. 2. Evolution of (3 at the inner edge of the disk for the 
four runs starting from (3 — 0.5 (dotted), 1.0 (solid), 5.0 
(dashed) and 10.0 (dot-dashed). Time is given in terms of 
the orbital period of the inner radius. 



Fig. 3. Comparison of runs with different azimuthal res- 
olutions and approximations of the Poisson kernel, for 
(3 — 1.0. The radial resolution is fixed at 256 grid zones. 
The dashed curve shows the results of the run with 128 
azimuthal grid zones and the standard kernel; dotted is 
for 32 zones and the smoothed kernel. 

Figure [l] shows the evolution of the rms velocity within 
the disk for values of (3 = 0.1, 0.5, 1.0, 5.0 and 10.0. As is 
shown, the perturbed velocity grows rapidly for (3 = 0.5 
and (3 = 1.0, exceeding the initial random noise after ap- 
proximately 10 orbits. The growth is approximately ex- 
ponential during this time. We will show in Sections i.l 
and 4.6 that the instability appears as a spiral mode, in 



agreement with the theory. For the other values of (3 we see 
very little action in comparison. The (3 — 0.1 and (3 — 10.0 
cases are virtually stable with the (3 = 5.0 being somewhat 
more active but failing to achieve the same amplitude as 
the most unstable cases. This run is still interesting to 
observe as the growth of the instability is slower and the 
evolution of the spiral wave is somewhat cleaner, as pre- 
sented later. 

Also shown in the figure is the rms magnetic field. This 
remains relatively constant over time for all cases, with 
only moderate increases in (3 occuring for the unstable 
cases {i.e. (3 = 0.5, 1.0 and 5.0) as shown in Figure |[ 
The maximum change in (3 occurs for the (3 — \.Q case. 
We see it decreasing to approximately (3 = 0.6 after the 
initial development of the instability as the field and mat- 
ter are accreted to the inner boundary. After this time 
we see some oscillatory features and finally a very gradual 
decrease to a value of (3 ~ 0.3 at around 200 orbits. Hence 
the plasma beta decreases to approximately 30% its initial 
value for the most active case. 

In the subsequent sections we examine more closely 
the most active values of plasma beta (namely (3 = 0.5, 
(3 = 1.0 and (3 — 5.0). As seen from the growth and de- 
velopment of the perturbed velocity, the cases (3 = 0.5 
and (3 = 1.0 are similar and this holds true for many 
other characteristics in the structural evolution. Similar 
features hold also for the (3 = 5.0 case with development 
being somewhat slower. 
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4-2. Effect of the spatial resolution 



As discussed in Section 3.3, the Poisson kernel given by 
Binney & Tremaine (1987), and used in the computation 
of the magnetic potential, is only valid for high azimuthal 
resolution (comparable with the radial one). For this rea- 
son the simulations presented in the majority of this paper 
are for 256 x 128 grid zones in u and 4> respectively. The 
high radial resolution is needed to obtain a good descrip- 
tion of the AEI, which depends critically on a resonance 
localised at the corotation of the modes (where their angu- 
lar phase velocity equals the rotation velocity of the gas). 
On the other hand the instability is predicted to occur 
only for very low azimuthal wavenumber m, so that a high 
angular resolution would be necessary only to describe in 
details the dissipation of wave energy at small scales, as 
the non-linear evolution leads to wave steepening an possi- 
bly to the formation of a spiral shock. This goes far beyond 
the purpose of the present paper and would require much 
more elaborate simulations. 

Thus simulations with low angular resolution, and in- 
creased radial resolution, would be enough for our present 



purposes. As discussed in Section 3.3, this leads to a nu- 



merical instability if we u se the regularisation, given by 
Binney & Tremaine ( 1987 ), of the logarithmic singularity 
of the Poisson kernel; we have made low azimuthal reso- 
lution runs with a different regularisation, corresponding 
to a smoothing by the finite vertical thickness of the disk. 
In this section we compare two simulations: one using the 
standard Binney & Tremaine kernel at 256 x 128 and an- 
other with a smoothed kernel (assuming a disk aspect ratio 
e = 0.1) with a resolution of 256 x 32. All other parameters 
are identical for the two runs, with f3 = 1.0. 

Figure |3| shows the different evolutions of the two cases. 
As seen, the initial evolution of the instability is virtually 
identical. A similar magnitude of u rms is achieved after 
approximated 20 orbits. Only the final, saturated state of 
the instability (which depends at least in part on the dis- 
sipation of wave energy at small scales) is slightly affected, 
both for the perturbed velocity and magnetic fields. This 
confirms that the smoothed kernel does give the correct 
physics we expect. 

4-3. Effect of the radial density gradient 

We now retain the same radial magnetic field profile, but 
change the density profile to check its effect on the growth 
of the instability. We have chosen not to try and repro- 
duce the profiles used in the linear analysis of TP99: these 
profiles, imposed by the numerical method used in TP99 
to solve the linear system, have a sharp and localised gra- 
dient at the corotation radius, in order to limit the range 
of density and/or magnetic field in the disk. They would 
thus be smeared rapidly by the non-linear evolution of the 
instability, precluding the more direct comparison of the- 
ory with numerical simulation one could expect. Work in 



Fig. 4. Comparison of runs with different density gradi- 
ents, for (3 = 1.0 at the inner edge of the disk. Lines 
plotted are dashed: E = const for a 256 x 128 grid, dot- 
ted: E oc r~ 3 / 2 for 256 x 128, dot-dashed: E oc r~ 3 / 2 for 
512 x 512 and dot-dot-dot-dashed: E oc r~2 for 512 x 512. 



progress with a different scheme for the solution of the lin- 
ear problem should relax this condition, and allow a direct 
comparison of the theory with the non-linear simulations. 

To study the stability criterion we chose two extra 
cases: one with a surface density varying as E ~ r -3/2 
and a steeper one with E ~ r~ 2 . These should be increas- 
ingly stable. The first case is identical to that chosen by 
Stehle & Spruit. According to the criterion (pi|), only the 
flat density profile should be subject to the AEI. 

Figure |i| displays the growth of the instability for these 
three different gradients of surface density, along with a 
run of different resolution. The initial random velocities 
vary slightly since the deviation from Keplerian rotation 
(thus the effect of the pressure gradient on the equilibrium 
rotation curve) is included in the rms values. 

This shows that, as expected, a radial surface density 
gradient stabilizes the modes. For E ~ r -2 the instability 
has totally disappeared. For E ~ r~ 3 / 2 it is still present, 
though weaker than with a flat density profile. Thus nu- 
merically the instability threshold seems to be steeper (in 
terms of the density profile) than given by the criterion 
(pij). This can be attributed to the fact that, as discussed 
by TP99, there are actually two different contributions to 
the growth rate: the corotation resonance gives the cri- 
terion (53). However a second mechanism (actually the 
first discovered by Tagger et at, 1990) corresponds, as for 



galactic spirals, to the emission of an outgoing spiral wave, 
emitted beyond the corotation radius. This outgoing spi- 
ral is very visible in the plots, figures 5-7. This mechanism 
does not depend on the radial profiles (we refer to TP99 
for a more complete discussion). Although we expect it 
to be weaker than the corotation resonance, it shifts the 
instability threshold leaving the E ~ r~ 3 / 2 case weakly 
unstable. The presence of the outgoing spiral will be dis- 
cussed later in more details. In order to verify the insta- 
bility of the E ~ r~ 3 / 2 profile, we have performed a very 
high (512 x 512) resolution run, giving results similar to 
the lower resolution one (though the growth rate is slightly 
lower, as shown in Table |l|). 

4-4- Growth rates 

The plots of rms velocity (Figures [l], [|and [I]), show that in 
general the instability grows exponentially until it reaches 
a saturated state. We give in Table [l] the growth rates in 
terms of f2o - the angular velocity at the inner radius - 
of all the runs presented here covering the different reso- 
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Fig. 5. Evolution of the spiral structure at regular intervals in time for the f3 = 5 case. Starting from an initial random 
state the instability structure evolves from an m — 3 (three-armed) spiral to an to — 2 and an to = 1. The plots show 
the radial velocity. Times shown are in rotation periods at the inner disk edge. 



lutions, (3 and radial density gradients. The growth rate 
indeed appears to be strongest for the flat density gradient 
with [3 — 1.0, and to almost totally disappear when [3 is 
one order of magnitude higher or lower. Similarly, as the 
density gradient gets steeper, the growth rate decreases 
and the case with £ ~ r~ 2 is totally stable. 



p 


Resolution 


a (E oc r~ a ) 


Growthrate /f2o 


0.1 


256 x 128 





0.0001 


0.5 


256 x 128 





0.0346 


1.0 


256 x 128 





0.0369 


5.0 


256 x 128 





0.0058 


10.0 


256 x 128 





0.0002 


1.0 


256 x 32 





0.0371 


1.0 


256 x 128 


3/2 


0.0173 


1.0 


512 x 512 


3/2 


0.0143 


1.0 


512 x 512 


2 


0.0000 



Table 1. Comparison of growth rates for all runs per- 
formed. 



These values are comparable with the theoretically 
predicted ones, as given in TP99. As explained in the pre- 
vious section, a detailed comparison is difficult, since TP99 
use quite unrealistic radial profiles, in order to limit the 
difficulties associated with their solution method. 

4-5. Evolution of the spiral wave 

The evolution of the instability into a well developed spiral 
wave is best seen for the (3 = 5.0 case. Figure || shows the 
evolution from the initial random velocity field at intervals 
of 20 orbits. The plots are of the radial velocity in the inner 
region of the disk. 

By t = 60 orbits (at the inner edge of the disk), an 
m = 3 (3-armed) spiral is clearly evident. This then de- 
velops over time to lower wavenumbers, eventually becom- 
ing a tightly wound m — 1 spiral after about 200 orbits. 
The disk goes through some intermediate stages for which 
a mixture of wavenumbers co-exist. This is discussed in 
more detail in the following section. 

The development of the spiral from higher to lower 
wavenumbers is again in agreement with the theory: the 
most unstable azimuthal wavenumber decreases from m = 
a few to to — 1 as (3 decreases. This decrease of f3, shown 
in Figure ||, occurs as magnetic flux is advected inward 
with the gas to the inner region of the disk. 

All the runs show similar behaviour to the (3 — 5.0 
case. The spiral always develops from high to small 
wavenumbers, eventually becoming to = 1 in all cases. 



Fig. 6. Contour plot of radial velocity as a function of ra- 
dius and time, at (f> = 0, in the inner region of the disk 
for the (3 — 1.0 case. Oblique features show a wave trav- 
elling outward, while horizontal ones indicate a standing 
pattern. 

Fig. 7. Outgoing waves propagating outward. Radial ve- 
locity is plotted radially at 4> = , in the (3 — 0.5 run. The 
times plotted are t = 100.0 (solid), t = 100.2 (dashed), 
t = 100.4 (dot-dashed), t = 100.6 (dot-dot-dot-dashed), 
t = 100.8 (dotted). 

Fig. 8. Vorticity in the inner disk region in the (3 = 0.5 
run. 

Fig. 9. Velocity field in the inner disk region, in the case 
shown in Figure |[ Large azimuthal velocities close to the 
inner edge have been removed from the plot for clarity. 

4-6. Time analysis of the development of the instability 

We further this analysis of the evolution of the instability 
by performing a spectral analysis of the development of 
its structure with time. 

Figure |^ shows a contour plot of the radial velocity, 
measured on the axis <j) = 0, as a function of radius and 
time, for the (3 — 1 run, with lighter regions indicating 
higher velocities. This provides some essential details on 
the nature of the instability. It becomes visible after t ~ 20 
orbits. Higher velocities are seen to occur close to the inner 
radius of the disk. In the innermost region (r 2), the con- 
tours shown are horizontal, whereas farther out they are 
oblique with a positive slope: the latter represent a wave 
travelling outward, while the former indicate a standing 
pattern. This is the structure expected for the AEI which 
(like self-gravity driven galactic spirals) distinguishes two 
regions in the disk, separated by the corotation radius 
(where the angular phase velocity of the wave equals the 
rotation frequency of the gas) : inward from corotation, the 
instability establishes a standing pattern formed of an out- 
going and an ingoing spiral. Beyond corotation, it emits a 
trailing spiral wave travelling outward. 

This all occurs very close to the centre of the disk. 
Virtually no effects arc observed for r much greater than 



Fig. 10. The same flow pattern as in Figure | in a 
Cartesian projection. The vortex is clearly visible around 
r ~ 1.5. 
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20. This is also a good indication that the outer boundary 
will have little effect on the results. Virtually no waves 
travel to the outer boundary and hence reflections and 
subsequent interactions can be neglected. Any significant 
signal returning from the outer boundary would appear in 
this plot as an oblique pattern with a negative slope. 

Figure shows in more details the outgoing charac- 
ter of the wave in the outer region. Figure || shows the 
vorticity in the flow. A sharp feature is seen at r ~ 1.5. 
This again is expected from the AEI, which grows by 
coupling the spiral pattern to a Rossby vortex it gener- 
ates at its corotation radius. The sharp vorticity feature 
is a manifestation of this vortex. This will play a crucial 
role in future studies: the instability mechanism lies in the 
fact that the spiral inward from corotation extracts energy 
and angular momentum from the flow (thus causing ac- 
cretion), and transfers them either to the Rossby vortex 
or to the outgoing spiral beyond corotation. The vortex 
is usually more efficient (as discussed in Section [D|). It 
can be considered as a twisting of the footpoints of mag- 
netic field lines in the disk. Thus, if we took into account 
a small but finite density above the disk, this twisting 
would propagate as Alfven waves, carrying vertically to- 
ward the corona (where they could power a jet or an out- 
flow) the energy and angular momentum extracted from 
the disk. This is the process which inspired the name of 
the Accretion-Ejection Instability. 

In order to show the vortex more clearly, we have plot- 
ted in Figure ^ the flow pattern in the inner region of the 
disk. A turning point in the flow is seen close to the inner 
boundary in the top-right quadrant of the plot. This vor- 
tex is more clearly shown in the flow pattern of Figure |l^. 
For clarity, both figures show the deviations from the un- 
derlying Keplerian flow. A flow pattern around a singular 
point is seen in both plots. This point is at the corotation 
radius. 

Figure [ll] shows details of the spatial structure of the 
waves during four intervals of time, for the (3 — 5.0 run, 
which shows most clearly the spectral evolution. For each 
time interval it shows two plots: the left one is a contour 
plot of the radial velocity, analogous to the one shown 
in Figure || The right plot shows contours of the Fourier 
transform of the left one, giving as a function of radius, 
the frequencies at which perturbations evolve during the 
given time interval. Frequencies of lu/Hq ~ 2.2, 1.5 and .7 
appear succesively, corresponding to the m = 3, 2 and 1 
spirals seen in Figure ||. These frequencies correspond to 
corotation radii which increase (uj/m decreases) as m de- 
creases, as expected from the theory. 

We note that results recently released by Hawley & 
Krolik ( 2000| ) and Hawley (2000b) show similar behaviour 
close to the inner edge of an accretion disk. Their numer- 
ical setup is much more complex than ours, and describes 
a 3D MHD model of the inner region of an accretion disk 
around a black hole with a pseudo-Newtonian potential. 
Although their goal is to study the magneto-rotational 



instability in a realistic 3D geometry, they note that af- 
ter some time a significant part of the radial transport 
is done by waves. They show, in the inner region of the 
disk, plots very similar to our Figure |^. A detailed com- 
parison is difficult, since this appears after an important 
evolution of the disk and magnetic field geometry. Fur- 
thermore, the relatively poor radial resolution imposed by 
the full 3D simulation would probably rule out a proper 
description of the AEI. However the similarity between 
their results and our Figure ^| (including the horizontal 
and oblique patterns discussed previously) is striking. We 
also note that, when this appears, the magnetic field has 
strongly evolved from their initial setup, as some vertical 
magnetic flux has been advected inward with the gas. In 
the inner region of the disk, they thus have an essentially 
vertical field of rather high amplitude (although the value 
of f3 is not given). We thus consider that the wave pat- 
terns they obtain could very well be a manifestation of 
the AEI, rather than waves generated non-linearly by the 
magneto-rotational turbulence. 

4-7. Accretion rates 

From the simulations we are able to estimate the efficiency 
of the instability in accreting matter towards the central 
object. At a given radius the accretion rate is calculated 
as: 



M(r) = 



2tt 



ru r (r, 0)£(r, <f)d<p. 



(26) 



However, this value is highly variable. Plotting against r 
at one point in time shows very little of the overall ef- 
fect of the instability. We therefore show values that are 
averaged over space or time, in order to obtain a better 
understanding of the overall accretion within the disk. All 
figures shown are for the (3 = 1.0 run. 

Figure |lj shows the accretion rate against time, av- 
eraged between r = 1 and r = 10. This region is chosen 
since the instability has little effect at larger radii. This 
shows that even though the accretion rate is highly in- 
termittent, it is predominantly positive. Three bursts of 
accretion are seen before the first 100 orbits, at t ~ 30, 
t ~ 50 and t ~ 75 orbits. These burst are also seen in the 
plots of rms velocity in Figure |l|. 

However, the physics of the instability lets us expect a 
more complex behavior: spiral density waves carry a neg- 
ative flux of energy and angular momentum inside their 
corotation radius, and a positive one beyond corotation: 
this means that a wave inside corotation grows by ex- 
tracting energy and angular momentum from the disk, 
and transfering them either to the Rossby vortex, at the 
corotation radius, or to the outgoing spiral, beyond it. 
Accordingly, we expect to see (at least during the initial 
phase when the instability grows following its linear be- 
havior) accretion in the inner part of the disk, and outward 
motion (negative M) beyond the corotation radius, which 
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Fig. 11. Contour plot of the radial velocity against time close in the inner disk region, in the f3 = 5.0 run, with 
associated spectral analysis. Time ranges shown are for 50 to 75 orbits, 100 to 125 orbits, 150 to 175 orbits and 200 
to 225 orbits. Spectra show first a feature at uj/tto ~ 2.2, then one at ~ 1.5, and finally at ~ .7. They correspond to 
the m = 3, 2 and 1 spirals succesively formed in the disk. 



Fig. 12. Accretion rate, averaged between r = 1 and r = 
10, against time for the (3 — 1.0 run. 



is typically of the order of 1.5 — 2. This is exactly what we 



see in Figure 13 1, which shows the accretion flux in the 



inner region of the disk, averaged between times 10 and 
25 and for the same run with (3 = 1. However the result is 
still quite noisy, especially at large radius: this turns out 
to come from our initialization, where we have introduced 
strong (~ .1) random fluctuations of the radial velocity, 
from which the instability grows. By the time these fluctu- 
ations have dissipated or traveled away, they have caused 
significant radial flows. In order to check this we have per- 
formed a new simulation with the initial rms noise reduced 
to .001: Figure [l3| b shows that the noise has strongly de- 
creased at large radius, while the inner region is qualita- 
tively similar (note that absolute values of M cannot be 
compared, since the instability now takes longer to grow 
from the initial conditions to a given amplitude). Finally 
we run a new case where the initial random perturbation 
now depends on r as a Gaussian centered around r = 5: 
Figure |l3| c shows the resultant accretion flow, which has 
the same structure as in Figure |l3|b, but where the fluc- 
tuations at large radius have totally disappeared. Figure 
|l3| d is a zoom on the inner disk region, showing that as ex- 
pected M is positive inside a radius r ~ 1.7, and negative 
beyond, going to zero at large radius. 

At larger times, however, accretion becomes positive 
beyond the corotation radius: this is shown by Figure |l4|, 
where M is averaged over the whole simulation. It is in- 
teresting to see that, in the region between corotation and 
r = 10, most of this positive accretion occurs by bursts. 
These (already noted in Figures [l] and [l2]) are shown on 
Figure [l5], where M is shown as a function of r and t. 
The bursts are seen as strong streaks propagating outward 
from the corotation region. We conclude that this is due to 
non-linear effects (e.g. shocks, or global re-arrangements 
of the current profile) not included in the linear theory, 
but responsible for the saturation of the instability after 
a finite time. We defer their full characterization to fu- 
ture work where more realistic magnetic field and current 
profiles will be used. 

It is customary to measure the efficiency of accretion 
in terms of a turbulent viscosity, parametrised with the 
parameter a ss of Shakura & Sunyaev (1973). This al- 



lows a comparison of the instability mechanism with stan- 
dard a-disk model. We thus convert the accretion rate to 
an equivalent value of a ss , although we emphasize that 
the physics of the AEI (which forms large-scale, quasi- 



stationary patterns where the accretion energy is carred 
away by waves) is very different from the assumptions un- 
derlying the Shakura-Sunyaev model, i.e. small-scale tur- 
bulence leading to local deposition of the accretion energy. 
Here, as explained in previous sections and in more de- 
tails in TP99, the instability mechanism relies essentially 
on the excitation of the Rossby vortex by the spiral wave. 
This occurs through the long-range action of the magnetic 
potential <&m- Thus the accretion energy is not deposited 
locally, but transported at large distances by the magnetic 
stresses figuring in equations (1-2). The present estimate 
can thus be taken only as a convenient measure of the 
efficiency of the instability to cause turbulent transport. 
Furthermore, a reliable assesment of the accretion rate due 
to the AEI would require taking into account the vertical 
emission of Alfven waves (i.e. include the effect of a finite 
density corona above the disk) , and effects associated with 
the finite thickness of the disk. 

From standard thin disk approximations of steady ac- 
cretion disks, one has 



uTi = — 



M 
37 
M 
3^' 



1 - 



1/2 



(27) 



(Frank, King & Raine 1992) where v is the viscosity, E the 
surface density and r* is the radius of the central object. 
The standard prescription for the viscosity is given by 



v = a ss c s h. 

Combining Equations ^ 
disk assumption, that c s 



s(r,t) 



M(r, t) 



Sneu^r, i)rE(r, t) 



(28) 

and [2^ along with another thin 
= hil, yields 

(29) 



where e is the disk aspect ratio. The values of u^(r, t) and 
E(r, t) represent azimuthally averaged quantities resulting 
in an a ss that depends on radius and time. 

Figure [l^ shows the effective a ss against time when 
averaged between r = 1 and r = 10, for the /3 = 1 run 
with Gaussian initial perturbations. It is negative while 
the wave is in its linear growth phase, but becomes positive 
in the non-linear stage. It is again strongly intermittent, 
showing the same bursts as in previous figures. 

It is useful to compare values of a ss with values de- 
rived from numerical models of turbulent accretion disks, 
for which the use of the a parameter is more applicable. 
The use of shearing-box models (e.g. Brandenburg et al. 
1995, Brandenburg et al. 1996a, Stone et al. 1996, Hawley 
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Fig. 13. The accretion rate M as a function of radius, for the (3 = 1 run, in the linear growth phase of the instability. 
Results are shown from runs with different initial conditions for the random velocity fluctuation in the disk, a (top 
left): strong and constant (constant meaning maximum Mach number), b (top right) weak and constant, c (lower left) 
weak, peaked around r — 5, and d (lower right) magnification of the inner region of c. As expected for spiral waves in 
the disk, accretion is positive inside the corotation radius, and negative outside. 



Fig. 14. The accretion rate M as a function of radius, 
for the [3 = 1 run, averaged over the whole simulation. 
Accretion has become positive up to r ~ 10. 

Fig. 15. Evolution of radial structure of accretion rate for 
(3 = 1 with initial perturbations given a Gaussian profile. 

Fig. 16. Estimation of an equivalent Shakura-Sunyaev a- 
parameter from the accretion rate, averaged between r = 1 
and 10, in the (3 = 1 run. 



et al. 1996) over recent years have given insights to the 
turbulent transport of energy via calculations of Reynolds 
and Maxwell stresses as well as other methods (such as 
from the overall increase in thermal energy) . Various val- 
ues of a calculated from different quantities by Branden- 
burg et al. (1996a) indicate that their a ss is typically of 



the order of 0.001 to 0.01 with maxima generally around 
0.01 to 0.02. These are comparable to those of Stone et al. 
( p9"6[ ) and and Hawley et al. ( p96| ). 

With the AEI simulation in 2d for (3 = 1.0 we see sim- 
ilar values of the viscosity parameter, a ss . The averages 
over time and space yield values typically between 0.001 
and 0.006. Maximum values at a given position and time 
are around 0.06. These are still somewhat smaller than 



the recent values given by Hawley & Krolik ( 200C ) for 
their 3d model of the inner accretion disk and also that of 



Armitage (|1998Q . 

The strength of the instability, depending on (3 and 
other parameters, affects the value of a ss . Table || shows 
the maximum values of a,, from each of the runs. 



(3 


Otss 


0.1 


0.011 


0.5 


0.122 


1.0 


0.056 


5.0 


0.030 


10.0 


0.013 



Table 2. Comparison of equivalent values of a sa derived 
from the accretion rate for each case of f3. 



Again we see that the instability is strongest for mod- 
erate values of (3, i.e. around j3 ~ 0.5. Indeed, we see 
that values of a are strongly peaked for a = 0.5 and are 
more comparable to values desirable for models of cata- 



clysmic variables (Cannizzo 1993 ). Recen t globa l simula- 
tions of the MRI (Armitage |l998| , Hawley |2000a|) indicate 
that values of a derived from global numerical simulations 
are generally higher tha n local shearing-box simu lations 
{e.g. Brandenburg et al. 1996a, Stone et al. 1996, Miller 



& Stone 200C) and hence the possibility exists of cata- 
clysmic variables occuring as a result of a combination of 
both the AEI and the MRI. We must also emphasize that 
the physics discussed here concerns only the innermost 
part of the disk. Presumably, if accretion proceeds faster 
at larger radii, the resulting inflow of matter and magnetic 
flux will strongly affect the evolution of the AEI and its 
consequences on inward transport. 

5. Conclusions 

In this paper we have presented a two-dimensional model 
of the inner region of an accretion disk threaded by an 
initially vertical magnetic field, surrounded by vacuum. 
Starting from an equilibrium situation with additional 
random small scale fluctuations, an instability is seen to 
occur on a timescale of a few tens of orbits. The instabil- 
ity appears as a low azimuthal wavenumber spiral wave, 
amplified mainly by (and feeding) a Rossby vortex. These 
results are comparable (with a slightly different physical 
setup and a more adapted numerical one) with the ones 
obtained recently by Stehle and Spruit ( 2000p . 

The instability is seen to occur most strongly for a 
plasma beta of f3 ~ 1 (i.e. equipartition between the gas 
thermal pressure and magnetic pressure in the disk) , with 
its amplitude and growth rate rapidly decreasing for both 
stronger and weaker fields. Also, the growth rates (and the 
existence of the instability) depend on the radial gradients 
of equilibrium quantities. 

As the spiral wave develops and affects the disk struc- 
ture (in particular decreasing its (3 value) it is seen to 
progress from larger to smaller azimuthal wavenumbers 
until it reaches m = 1, i.e. a one-armed spiral. 

These properties (timescales, development, character- 
istics, stability criterion) allow us to identify the insta- 
bility as the Accretion-Ejection Instability, and globally 
to confirm qualitatively and quantitatively the behaviours 
predicted theoretically by TP99. 

The instability is seen to generate accretion in the in- 
ner region of the disk. The amplitude of the accretion - 
calculated through an equivalent value of the standard 
a pa rameter of viscous transport (Frank, King & Raine 
1992 ) - is comparable, when compared to other numeri- 
cal simulations of accretion disks (e.g. Brandenburg et al. 
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1995| , Stone et al. 1996), with the one generated at lower 
magnetic pressure by the magneto-rotational instability 
(MRI) of Balbus & Hawley ( |l99l| ) . 

Furthermore, we have noted that waves observed in 
more recent publications of global, 3D simulations of the 



inner region of accretion disks (Hawley & Krolik 2000, 
Hawley 2000b| ) show properties (development of spiral 
waves with a characteristic propagation pattern) which 
are very similar to the ones described here. These waves, 
which seem to cause a significant part of the accretion, 
develop once the global magnetic structure has evolved 
in the inner region of the disk to become quite similar to 
the one used here. It is thus quite possible that the AEI 
is also present in these simulations, once the Magneto- 
Rotational Instability has caused the disk to evolve to a 
favorable magnetic configuration (in terms of geometry 
and strength) in its inner region. 

Future work needs to extend the current simulations 
to 3D, in two steps: the first one will be to consider a disk 
which for simplicity is still infinitely thin, but embedded 
in a corona of small but finite density: this will allow us to 
describe the vertical ejection, as Alfven waves, of the en- 
ergy and angular momentum extracted from the disk and 
presently stored in a Rossby vortex. It will thus be pos- 
sible to consider how they can ultimately be deposited in 
the corona, where they could energize a jet or an outflow. 

A second step will be to turn to fully 3D simulations. 
This leads to numerical constraints which are quite dif- 
ferent from the ones encountered in the simulation of 
magneto-rotational turbulence, since the physics of the 
AEI is quite different {e.g. it needs very few grid points 
within the disk in the vertical direction, but a large num- 
ber in the radial one). This will allow us to address the 
acceleration of the gas, across a slow magnetosonic point, 
to finally form a jet. 
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